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ABSTRACT 

Recent focus on the importance of cold, unshocked gas accretion in galaxy formation - 
not explicitly included in semi-analytic studies - motivates the following detailed com- 
parison between two inherently different modelling techniques: direct hydrodynamical 
simulation and semi-analytic modelling. By analysing the physical assumptions built 
into the gasoline simulation, formulae for the emergent behaviour are derived which 
allow immediate and accurate translation of these assumptions to the GALFORM semi- 
analytic model. The simulated halo merger history is then extracted and evolved using 
these equivalent equations, predicting a strikingly similar galactic system. This exer- 
cise demonstrates that it is the initial conditions and physical assumptions which are 
responsible for the predicted evolution, not the choice of modelling technique. On this 
level playing field, a previously published GALFORM model is applied (including addi- 
tional physics such as chemical enrichment and feedback from active galactic nuclei) 
which leads to starkly different predictions. 



1 INTRODUCTION 



Two different approaches are widely used to test theories 
of galaxy formation. Both make use of developing compu- 
tational resources to integrate a set of coupled differential 
equations forward in time, where each equation applies phys- 
ical constraints, or empirical laws, to some of the system's 
properties. 

The semi-analytic method (White & Frenk 1991) can 
be characterised by the intent to encapsulate as much of the 
physical behaviour as possible in the equations themselves, 
chosen to describe systematic properties (virialisation, con- 
servation of angular momentum, etc.) and thus leaving the 
minimum to numerical calculation. 

Inherently different to this is the approach taken by 
simulations, where the ethos is to apply equations which are 
known to govern the particular properties (e.g. the gravi- 
tational force between two particles) and demonstrate, by 
skilfully programmed numerical calculation, that these re- 
produce macroscopic properties. 

In practice, adjustment will always need to be made to 
account for the fact that each simulated particle represents 
many, many real particles. In traditional simulations which 
considered only gravity, this adjustment is purely mathemat- 
ical. When gravity is the dominant interaction, results are 
therefore very robust and simulations are mostly accepted 
as being a faithful realisation of the Universe at large scales, 
where this is the case. A notable measure of this acceptance 
is that semi-analytic models have taken to using samples 
of halos taken directly from simulations (Kauffmann et al. 
1999; Helly et al. 2003; Hatton et al. 2003), rather than sam- 
ples created statistically (Bower et al. 1993; Lacey & Cole 



1993) using methods based on the analytic arguments of 
Press & Schechter (1974). 

In general, though, there must come a scale below which 
the governing equations are no longer describing particle in- 
teractions but emergent phenomena. At these smaller scales, 
the two largely contrasting approaches begin to look more 
similar. As an immediate example, simulations which en- 
compass an entire galaxy have not attempted to resolve in- 
dividual stars and must therefore apply an equation which 
approximates the conversion rate of large bodies of gas into 
stars, based on locally averaged density. An equivalent and 
sometimes identical equation appears in semi-analytic mod- 



els (1 3.71 



Both approaches have had some success in accounting 
for the observed properties of galaxies. Predictive power, on 
the other hand, relies on surety in the underlying equations, 
or a clear understanding of how differing underlying assump- 
tions map on to final observable quantities. 

To bolster this crucial understanding, it is extremely 
helpful to study, in detail, the predicted evolution of a single, 
large galaxy and its satellites, as followed by these two differ- 
ent techniques. If the analytic model begins from the precise 
halo merger history found in the simulation, and applies the 
same assumptions about the physics at smaller scales, will 
it predict the same observable properties and, if not, where 
and why do the two calculations diverge? 



1.1 Previous research 

Existing quantitative comparisons between simulations and 
semi-analytic models have focussed on samples of galaxies. 
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The approach which is adopted here - adapting the semi 
analytic model to emulate the assumptions made by the 
simulation — was also used by Helly et al. (2003) to demon- 
strate that the cooled mass predicted by the two approaches 
can be consistent. In particular, consistency was demon- 
strated over two orders of magnitude in halo mass. Earlier 
works by Benson et al. (2001) and Yoshida et al. (2002), us- 
ing similar techniques, reached broadly similar conclusions. 

A broader comparison was made by Cattaneo et al. 
(2007) of two predicted galaxy populations: one by an SPH 
model and the other by version of the GallCS semi-analytic 
model, modified to emulate the SPH method. Again, consis- 
tency is demonstrated between the two techniques once the 
semi-analytic model is appropriately "stripped down", but 
this consistency comes at the expense of agreement with 
observational constraints, which has been met by the full 
version of GallCS. 



Table 1. Resolution details for the gasoline simulation 



1.2 Outline 

To extend this previous research, in this present work we 
will compare in detail the formation of a single galaxy in 
a gasoline numerical simulation and in the galform 
semi-analytic model. The structure of this paper is as 
follows: 

§2: A review of the broader aspects of the gasoline code 
and of this particular simulation. 

§3: Detailed analysis of the way that gasoline models each 
of the following key physical processes, together with an ex- 
planation of how the particular approach in question can 
be emulated in the galform model. 

§3.1: Merging 

§3.2: Accretion 

§3.3: UV heating 

§3.4: Gas distribution 

§3.5: Cooling 

§3.6: Disk formation 

§3.7: Star formation 

§3.8: Feedback 

§3.9: Disk stability 
Particularly relevant aspects of the galform model are 
reviewed (or referenced) when appropriate. General aspects 
of the model can be found in Cole et al. (2000) and Bower 
et al. (2006). 

§4.1: By applying these same underlying physical assump- 
tions to galform, the two modelling techniques are shown 
to make broadly consistent predictions. 

§4.2: Having demonstrated this consistency, galform is 
then used to find how this halo would evolve under the 
physical assumptions made in previous studies (which have 
been shown to successfully match the collective properties 
of galaxies). 



Dark Matter particles: 


7.t 


x 10 5 M Q 


Gas particles: 




x 10 5 Af Q 


Star particles: 


3.8 


x 10 4 Af o 


Force resolution: 


0.2 


kpc 



2 GASOLINE 

The "MW1" simulation analyzed in this paper was gener- 
ated using Gasoline, an N-Body+Smoothed Particle Hy- 
drodynamics (SPH) Parallel Treecode (Wadsley, Stadel & 
Quinn 2004), within a flat, A-dominated WMAP1 cosmology 
(Qo = 0.3, A=0.7, h = 0.7, a 8 = 0.9, Q b = 0.039). A mod- 
erate resolution version of this run was originally discussed 
in detail in Governato et al. (2007), and the higher resolu- 
tion version analised here has subsequently been used for a 
number of studies (Zolotov et al. 2009; Read et al. 2009). 
High mass and force resolution (see Table [T| was achieved 
using the volume renormalization technique (Katz & White 
1993). 

This simulated galaxy matches the observed mass- 
metallicity relationship at varying redshifts (Brooks et al. 
2007; Maiolino et al. 2008). This success demonstrates that 
these simulations overcome the historic "overcooling prob- 
lem" (Mayer, Governato & Kaufmann 2008), in which gas 
in simulations cools rapidly, forming stars too quickly and 
early, and thus producing a mass-met allicity relation that 
is too enriched at a given stellar mass. Brooks et al. (2007) 
showed that reproducing the mass-metallicity relationship 
was a result of the adopted feedback prescription. This feed- 
back mechanism, combined with high numerical resolutiorj^] 
avoids spurious loss of angular momentum. As a conse- 
quence, simulated galaxies also match the observed Tully- 
Fisher relationship (Mayer, Governato & Kaufmann 2008; 
Governato, Mayer & Brook 2008; Governato et al. 2009). 

Gasoline's star formation and supernovae feedback 
scheme is described in detail in Stinson et al. (2006), and 
parameters from Governato et al. (2007) are adopted for 
MW1. The star formation recipe reproduces the Kennicutt- 
Schmidt law, and each star particle is born with a Kroupa 
initial mass function (Kroupa 1993). Combined with a pre- 
scription for a cosmic UV background (Haardt & Madau 
1996) that mimics reionization and unbinds non-collapsed 
baryons from simulated halos with total masses below about 
10 9 Mq, this feedback scheme drastically reduces the number 
of luminous satellites containing a significant mass of stars, 
avoiding the well known "substructure problem" (White & 
Frenk 1991; Kauffmann, White & Guiderdoni 1993; Quinn, 
Katz & Efstathiou 1996; Moore et al. 1999). 

As described in Stinson et al. (2006), supernova energy 
is deposited into the ISM mimicking the blast-wave phase of 
a supernova (following the Sedov- Taylor solution). Cooling 
is turned off for gas particles within the supernova blastwave 
radius for the duration that the remnant is expected to adi- 



abatically expand (£3.8 1 . Supernova feedback regulates star 
formation efficiency as a function of halo mass, resulting 
in the mass-metallicity relationship described above. This 
regulation of star formation also leads to realistic trends in 



§5: Summary. 



With ~ 4.8 X 10 total particles within the virial radius at z=0. 
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gas fractions, with the lowest mass galaxies being the most 
gas rich (Brooks et al. 2007), and reproducing the observed 
incidence rate and column densities of Dampled Lyman a 
systems at z = 3 (Pontzen et al. 2008). 

MW1 and its most massive progenitor halo were identi- 
fied at each simulation output step using AHF0(Gill, Knebe 
& Gibson 2004; Knollmann & Knebe 2009). AHF determines 
the virial radius, r v , for each halo at each output time step 
using the overdensity criterion for a flat universe following 
Gross (1997). The full gas accretion history of MW1 is de- 
scribed in detail in Brooks et al. (2009), and used in this 
study for comparison to Galform. 



3 TRANSLATION FROM SIMULATION TO 
SEMI- ANALYTICS 

3.1 The merger tree 

The common point between the two different models of 
galaxy formation is the halo merger tree. This is gener- 
ated by finding bound structures in the "MW1" simulation 
of Brooks et al. (2009) using AHF. A virial mass is assigned 
to each structure and its fate at each timestep is categorised, 
being either: 

1. Isolated; 

2. Inside the virial radius of a larger system; or, 

3. Merged. 

Two criteria are used to classify a galaxy as "merged": 

(a) If the DM mass decreases by more than 50% in a single 
timestep; or, 

(b) If the number of DM particles falls below 64. 

These are both significantly different from the definition 
adopted by galform (see Cole et al. 2000), so a detailed 
comparison of merger times, along the lines of (Jiang et al. 
2008), is deferred for another study. 

When a structure changes from category 1 to category 
2, the galform code updates its status to a "satellite" and 
predicts its subsequent evolution based on its approach ve- 
locity, v and position r at this juncture. These two proper- 



ties are combined into the dimensionless parameter (i Al I, 



which characterises the orbit. The distribution of values for 
O in this simulation was found to be consistent with the 
most recent assumption applied in galform (Fig. [TJ . 

This orbital information is then used to estimate a 
merger time based on the standard Chandrasekhar formula 
( |A6| | as defined by Lacey & Cole (1993). Details of this part 
of the model can be found in Appendix [A] and in Cole et al. ( 
2000; §4.3). 

The resulting dark and baryonic components are shown 
in Fig. [2] Being generated purely from the derived merger 
tree, the two realisations might be expected to be identical. 
One reason this is not quite the case is because the total 
virial mass, M v , is measured directly, but the relative com- 
position of baryons and dark matter may vary a little from 
halo to halo, and from one realisation to the other. 



2 AMIGA'S Halo Finder, available for download at 
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Cole et al. (2000)_ 
GASOLINE 




Figure 1. Orbital parameters randomly selected from the distri- 
bution (A5) of Cole et al. (2000) are shown in the dashed line, 
while the dotted line shows the distribution from Benson (2005). 
The solid line shows the distribution found in the simulation of 
Brooks et al. (2009) (which results from applying equation ( |A4| 
to the positions and velocities of the satellites). 



Previously, galform required that the halo mass was 
monotonically increasing as a function of time along each 
branch of the merger tree (i.e. each halo should be more 
massive than the sum of the masses of its progenitor halos). 
While this condition is always fulfilled for merger trees gen- 
erated using the extended Press-Schechter theory, it is not 
necessarily fulfilled for merger trees extracted from N-body 
simulations. 

For this reason, previous use of N-body merger trees in 
galform has required some adjustment of halo masses^] to 
enforce this monotonic behaviour (Helly et al. 2003). Un- 
fortunately, this adjustment is somewhat arbitrary and can 
lead to significant changes in halo mass when adjustments 
are propagated through the merger tree. 

Monotonicity in halo masses along branches of the 
merger tree is not enforced in this work. Instead, they are 
allowed to follow whatever evolution the N-body simula- 
tion provides. This requires some modification of the way 
in which galform assigns baryonic masses to halos. Previ- 
ously, each halo was assigned an initial mass of hot gas equal 
to: 



M h 



fib 



M v (new) — S M v (progenitor) 



progenitors 



(1) 



|http : //popia . ft ■ uam ■ es/ AMIGA 



where M v was the total mass of the halo and "progenitors" 
refers to halos that exist on a previous timestep but have, 
at some point, merged into the current "parent" (or its most 
massive progenitor). If halo masses do not monotonically 
increase, this could lead to negative gas masses. To ensure 
against this, the initial mass of gas is assigned to be the 
greater of and zero. 



3 The standard approach which was taken begins at the earliest 
times in the tree, checks whether each halo is more massive than 
the sum of the masses of its progenitor halos and, if it is not, sets 
its mass to be very slightly larger than that sum. This process 
was repeated for all halos in the tree. 
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Figure 2. Confirmation that the mass evolution of the "MW1" 
system which emerged in the gasoline simulation of Brooks et 
al. (2009) (solid lines) matches that contained in the extracted 
merger history used by galporm (dashed lines). 



The results are identical to the original implementation 
in cases where the halo masses are monotonically increasing, 
but this can result in the total baryonic mass not being pre- 
cisely equal to the fraction (Q^/Hm) of the total mass of the 
final halo. Fortunately, this difference is found to be small 
and, in any case, the simulated halos do not have precisely 
this ratio either. 



3.2 Accretion 

The temperature history of accreted gas was investegated by 
Keres et al. (2005) in a simulation containing 1120 galaxies. 
The distribution of gas particles in terms of their maximum 
past temperature is bimodal, suggesting that classification 
into hot and cold modes of accretion is appropriate. The 
contribution from each of these modes is then found to de- 
pend strongly on halo mass and on environment. Notably, 
the halo mass (M v ~ 3 x 10 n MQ) which separates low- 
mass galaxies dominated by cold accretion and higher mass 
galaxies dominated by the hot mode, is close to that found 
by analytic arguments (Birnboim & Dekel 2003). 

The decomposition of the accreting gas into different 
modes is a principle feature of the analysis by Brooks et al. 
(2009). The three modes of accretion are defined as follows 
(along with their analogues in the semi-analytic realisation): 

(a) Clumpy acccretion is defined to include "any parti- 
cles that have, at any previous output step, belonged to a 
galaxy halo other than the main halo we are considering." 
(Brooks et al. 2009) 

The information needed to evaulate this component in 
the semi-analytic model is readily available from the merger 
tree, being the virial masses, M v , of the existing, progenitor 
halos: 



M c i umpy (new) = 



fib 



progenitors 



AI V (progenitor) 



(2) 



The subsequent fate of any baryons within these substruc- 
tures is decided by their inbound trajectory and consider- 
ations of dynamical friction, as described in §3.1| and Ap- 
pendix [X] 

(b) Smooth accretion: Shocked component 




Time / Gyr 

Figure 3. An analysis of accretion onto the MW1 galaxy in gaso- 
line, in the light of assumptions associated with galporm. The 
top panel shows the cooling time for a gas particle at the virial ra- 
dius (JsJ at two relevant temperatures, together with the timescale 
for freefall under gravity J6J|. In the lower panel, the solid lines 
show how accretion is divided amongst these three modes in the 
simulation. The shaded areas show how the equivalent compo- 
nents would be estimated in galform. As a simple attempt to 
emulate |3^,), shocked and unshocked are distinguished by com- 
paring the cooling time at T = 9 /iqT v to the the freefall time 
(as explained in |3.2b). 



The following criteria are required for a gas particle in the 
simulation to be categorised as shocked: 



(a) T > -Tv 



and (6) Ap > 4 



(3%y 



V 3T J 

where 7v is the virial temperature defined in gasoline 

GM v /l)TlH 



3r v 



(3) 



(4) 



and Ap is the fractional increase in density since the pre- 
vious timestep. The latter is an entropy criterion and the 
factor of 4 is motivated by the Rankine-Huginot shock jump 
conditions for a singular isothermal sphere. Though low res- 
olution can artificially broaden shocks, Brooks et al. (2009) 
verify that this run is of sufficiently high resolution to resolve 
shock discontinuities. 

The virial temperature in galform is denned differently, 
being a factor of 3/2 higher, T v = |7^, so |3k) corresponds 
to T > ^T v . 
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Finding an appropriate analogue to this condition within 
the framework of galform is difficult. However, one relevant 
calulation that is made compares the cooling timescale for 
gas at temperature T: 

2p gas (r)A(T, Z) ' 

and the free fall time from a radius, r, enclosing mean den- 
sity p: 



tcoolM = 



(5) 



8nGp(r) 



(6) 



The cooling function A is found for given temperature T 
and metallicity, Z, by reference to the published tables of 
Sutherland & Dopita (1993). /9 gas is the gas density, p(r) 
the mean density of all the mass enclosed by radius r, and 
/imH is the mean particle mass. 

So, in the context of galform accreting gas is designated 
as shocked if: 

tcooi(7v, ?v) > t s (%,r v ) , (7) 
this being a reasonable, simple equivalent to |3^). 

(c) Smooth accretion: Unshocked component 

Gas in the simulation which fails to meet either of 
the criteria in |3| is categorised as unshocked. Similarly, 
accreting gas in galform is designated as unshocked if 

tcool < tff. 

The relative contributions of these three components 
are shown in Fig. [3] for both realisations. Being determined 
directly from the merger tree, it is unsurprising that the 
clumpy components are quite consistent. 

That the shocked component is overestimated by the 
semi-analytic model is entirely as anticipated by Brooks et 
al. (2009), though the eventual impact this has on the pre- 
dicted properties of the galaxy may not be quite as severe, 
as will be seen in the forthcoming sections. 

3.2.1 Alternative criterion 

Birnboim & Dekel (2003) investigate the conditions for the 
existence of virial shocks in galactic halos, culminating in 
the following expression for the existence of a shock: 



r v A(T v ) 



4 (pm H ) 2 \u \ 3 



< 0.0126 



2 

«0 = 



16fc B T v 
3/iniH 



(8) 



Referring to equations Q, (|5| and |6| above, this condition 
can be re-written: 

ttf < 0.2 t cool . (9) 

As can be seen from Figj3] the cooling time does not exceed 
the free-fall time by a factor of 5 until about 8 Gyr. The cri- 
terion |9| would therefore imply that all smoothly accreted 
gas will be unshocked before that point, corresponding to 
an unshocked:shocked ratio of 2 or 3, as opposed to 1/2 in 
gasoline or 1/6 according to Q. 

So, for this particular mass aggregation history, one sim- 
ple estimate Q, overestimates the shocked component in the 
simulation and another, more thorough estimate gives 
an underestimate. This suggests that the error in assuming 



spherical symmetry may unfortunately exceed any accuracy 
that can be gained by analysis of shock conditions. Calibrat- 
ing some of these analytic models to simulations may help 
quantitatively account for the complex geometry. 

3.3 UV heating 

A UV background is incorporated by Brooks et al. (2009) 
which is based on an updated model by Haardt & Madau 
(1996). This background "turns on at z — 6" and heats the 
gas in halos so that only those with " 100 or more dark 
matter particles are of sufficient mass to retain bound gas 
particles". For the resolution in this case, 100 DM particles 
corresponds to M cut = 7.6 x 10 7 M Q . 

The galform model is build with two corresponding 
parameters: A redshift, z cut , below which the UV back- 
ground heating is effective (set to z cut = 6 for consistency 
with the above) and a halo velocity, v cut below which all 
cooling is suppressed. For a virial overdensity, A v = 300, 
the mass cut-off, M cut , mentioned above translates to a cut- 
off halo velocity of: 

v cut = (A v /2)5 (1 + z)i (O M GM cut H )5 » 30 km s" 1 (10) 

These values are duly applied, but one key difference be- 
tween the modelling of this process remains. The UV back- 
ground in the simulation can slow cooling in larger halos 
without necessarily completely suppressing it, while in gal- 
form suppression is all or nothing (below and above v cu t 
respectively) . 



3.4 Hot gas distribution 

The distribution of hot gas that is assumed in the galform 
model was recently updated, in the light of work by Sharma 
& Steinmetz (2005), to the following density profile for hot 
gas of mass Mh: 



p(r) = 



P(0) 



(l + r/r H r 



P(0) 



M H /47rr H 



(11) 



This form was assumed in recent work on disk formation 
and structure (Stringer & Benson 2007), and will also be 
used here. 

The reason for this choice is clear from Fig. [4] which 
demonstrates that it is a satisfactory description of the dis- 
tribution of diffuse gas in the outer parts of this simulated 
system. In this Figure, the analytic form for the hot halo gas 
( |11[ ) and cold disk gas ( |21[ ) have been fitted to the simulation 
results, yielding the respective scalelengths: 



• Halo core radius, ru- 

• Disk scalelength, td- 



(with the initial constraint on each line that it integrates 
to the correct total mass for that component.) These scale- 
lengths, found from fitting the simulation results, can be 
scrutinised against those deduced from the usual analytic 
assumptions. 

Firstly, the halo core radius is typically assumed to be 
a constant fraction of the virial radius, which is not quite 
the case in this particular simulation (Fig. Q. However, in 
the absence of any obvious pattern in the evolution of th, 
the usual simple scaling is applied, with the choice: 



6 M. J. Stringer, A. M. Brooks, A. J. Benson, F. Governato 




Figure 4. The total gas density profile from the simulation of Brooks et al. (2009) at six different redshifts, shown as solid lines. The 
dotted lines show the analytic forms which are assumed by the galform model for the disk gas (eqn. |21| blue line) and diffuse halo gas 
(eq. [IT] red line). The sum total of these two components is also shown (magenta line). Each of these depend on a single parameter, 
fitted respectively by the disk scalelength, td and the core radius of the hot gas profile, rjj. (The integral under each analytic component 
is already fixed to match the total hot and cold gas mass within the virial radius.) Because there are so many more points to fit the hot 
component, these are much better matches than the cold gas. As it is the hot scalelength (not the cold gas scalelength) which is actually 
an input to galform, this convenient bias has not been corrected. 



r H = r v /3 , (12) 

being the best approximation. Fig. [5] compares this to val- 
ues found from fitting to the true gas distributions in Fig. 
[4] the result looking favourable. This straightforward scal- 
ing with the virial radius ( |12[ ) will therefore be used in the 
matched galform realisation. The cold gas scalelengths, 



also shown in Fig[5] are not inputs to galform, but are 
calculated within the model as described in §3.6| 



To understand the relevance of the hot gas distribution 
on the system's evolution as a whole, it is worth reviewing 
the precise treatment of gas cooling in galform. 
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Figure 5. The evolution of three key scalelenghts of the sys- 
tem: Virial radius, r v , halo gas scalelength, rjj and disk 
scalelength, r-Q l |21| . Dashed lines show the values that will be 
assumed by galform. Points show the values which best fit the 
gas profile of the simulation (Fig.[4]|. 



3.5 Cooling 

The cooling model described by Cole et al. (2000) determines 
the mass of gas able to cool in any timestep by following 
the propagation of the cooling radius in a notional hot gas 
density profil^] This profile is fixed when a halo is flagged 
as "forming" and is only updated when the halo undergoes 
another formation event. The mass of gas able to cool in any 
given timestep is equal to the mass of gas in this notional 
profile between the cooling radius at the present step and 
that at the previous step. The cooling time is assumed to 
be the time since the formation event of the halo. Any gas 
which is reheated into or accreted by the halo is ignored 
until the next formation event, at which point it is added to 
the hot gas profile of the newly formed halo. The notional 
profile is constructed using the properties (e.g. scale radius, 
virial temperature etc.) of the halo at the formation event 
and retains a fixed metallicity throughout, corresponding to 
the metallicity of the hot gas in the halo at the formation 
event. 

This work makes use of a new cooling model, which 
will be described in full detail in Benson et al. (in prep.). 
Rather than arbitrary "formation" events, this model uses 
a continuously updating estimate of cooling time and halo 



properties. The properties described in { 3.4 (density normal- 
ization, core radius) are reset at each timestep. The previous 
infall radiuaj^] (i.e. the radius within which gas was allowed 
to infall and accrete onto the galaxy) is computed by finding 
the radius which encloses the mass previously removed from 
the hot component in the current notional profile. 

The new model must supply an alternative estimate 

4 We refer to this as a "notional" profile since it is taken to repre- 
sent the profile before any cooling can occur. Once some cooling 
occurs presumably the actual profile adjusts in some way to re- 
spond to this and so will no longer look like the notional profile, 
even outside of the cooling radius. 

5 defined as the lesser of the cooling and freefall radii. 



of the time available for cooling in the halo, t ava ii, from 
which the cooling radius can be computed in the usual way 
(i.e. by finding the radius in the notional profile at which 
^cooi = tavail)- This is done by considering the energy radia- 
tion rate per particle, e r and the thermal energy per particle, 
ethj which are estimated by making standard assumptions: 



£th = -k B T v 



and 



e r =A(T v ,Z)n. (13) 

In terms of these quantities, the cooling time |5| is simply: 
eth(t) 



toooi(r,t) = 



£r (^"cool j t) 



(14) 



At any time, the model needs to identify some radius in the 
hot halo, r coo i, where the gas has had just enough time to 
radiate all its thermal energy. This as yet undetermined ra- 
dius is defined to satisfy the condition that the cooling time 
for gas particles at this point is equal to the time available 
for them to cool, which is estimated in terms of the mean 
energy radiation rate per partickrj e r : 



^avail(^) — 



J* e T (ti)dti 



(16) 



Equating (14 1 and (16) then gives: 

^cool — ^ avail 

§fc B T v _ /„' A(ti)n{ti)dti 



A(t)n(r cooh t) 



A(t)n(t) 



(17) 



where n is the mean number density inside the virial radius. 

This approach must account for the hierarchical assem- 
bly of a halo, so the denominator in the right hand side 
strictly includes a sum over all progenitors, i, of this system 
which exist at time ti. Incorporating this sum into (171, and 
rearranging, then yields a final expression to be solved for 
the cut-off density, p coo i = p(r coo i): 



|fe B T v 



(18) 



Pcooi (t) _ 

~ Jo [Ei A *(*0"i(*i)] dfc 

This equation is notably independent of halo structure, 
a factor which only comes in when the limiting density is 
used to compute the total cooled mass. The scalelength, rH, 



appears in the normalisation of the hot gas profile (111, and 
hence effects the mass enclosed by the cut-off density, p coo i. 



This can be appreciated graphically. Solutions to ( 18 1 
are illustrated in Fig. [6] for three choices of rn and three 
examples of cut-off density. For the higher cut-off density, 
a smaller rH (more centrally concentrated) will lead to a 
greater cooled mass. For the case of a low cut-off density, a 
centrally concentrated profile will actually lead to a lesser 
cooled mass, contrary, perhaps, t 



6 This is proportional to the mean number density n: 



N 



■ o(rH)An 



(15) 



The factor a does depend on halo structure (i.e. on rn) but, if the 
distribution of hot gas relative to the virial radius is the same for 
each halo, as implied in this case by \12\ , then this factor cancels 
out in JT7l. 
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Figure 6. Three normalised hot gas distributions (solid lines), 
corresponding to three different values of the scalelength rjr ap- 
plied to the profile Intersections between a solid line and 
one of the dotted lines correspond to solutions to equation ( |18| 
for a particular value of rjr and a particular cut-off density. The 
cooled fraction, M coo i /Afhot > is given by the integral under the 
solid curve up to the intersection with the relevant cut-off density. 



3.6 Disk formation 

In galform, galactic disks are assumed to form from the 
cooled gas, conserving the net angular momentum that it 
possessed when it was part of the halo. The assumed dis- 
tribution for the specific angular momentum, j, of the hot 
halo gas is taken from Sharma & Steinmetz (2005): 



dM(j ) 



oc j 



exp 



(19) 



The value of a was found to lie between 0.5 and 1.5, but the 
mean value a = 0.89 is adopted here. This is normalised to 
enclose the correct total mass and the mean specific angular 
momentum of the hot gas, (j), which is found directly from 
the spin parameter, 



A = 



(j) \E 



1/2 



GM 



3/2 



0') = 



1 



jdM(j) 



(20) 



Details of the calculation of the halo energy, E, can be found 
in Cole et al. ^000; §3.2.1). The disk scalelength, td, is re- 
calculated after each timestep such that angular momentum 
is conserved in the system. The calculation is based on the 
following surface density distribution (applied to each disk 
component, i): 



Ei{R) = Ei(0)e ''D 



[so Mi = 27rEi(0)r£] (21) 



and a circular speed, v c , which is constant with radius. The 
subsequent evolution of the disk is then followed by apply- 
ing physical assumptions which are as close as possible to 
those made in the simulation, as is explained in the following 
sections. 

As with the sub-halo trajectories, the parameter A can 
ordinarily be assigned at random, from a specified distri- 
bution. In this case, its value is measured for each of the 
simulated halos and passed on to galform, along with the 
rest of the merger tree information (£3.1 1. 

The goal here is not to prove that parameters like A, and 
O must necesarily be set by simulations in order to validate 
the results of any semi-analytic model. Rather, it is an ex- 
ample of the correspondence between these initial conditions 



and the properties of the galaxy which they produce. This 
should reassure theorists that, when such models are applied 
to large samples of galaxies, their individual properties will 
be representative as long as these governing, random, but 
physical parameters are drawn from the correct natural dis- 
tributions. 



3.7 Star formation 

At face value, gasoline and galform would appear to be 
applying the same assumptions regarding the conversion of 
cold gas to stars. These respective assumptions can be found 
from equation (4) of Stinson et al. (2006) and equations (4.4) 
and (4.14) of Cole et al. (2000Q 

dp* _ Pcoid 

dt tdyn 



(a) 



and 



(6) V = e* A ^ -(22) 

Tdisk 



where ip is the total instantanoues star formation rate in the 
whole disk, and p* is the local value per unit volume. The 
timescales which appear in these expressions have similar 
physical motivation, but are not identical. One, idyn, is the 
time for gravitational collapse of a spherically symmetric 
region of local density p co id- The other, Tdisk, is the angular 
period of the whole system at the half-mass radius, r i : 



(a) t d y 



\/ Gpcoh 



and 



(b) 



Tdisk 



n/2 

V c 



(23) 



So, despite the apparent similarity of the two expressions in 
( 22 1, they can lead to quite different instantaneous total star 
formation rates, ip. To illustrate this, one can consider the 
surface density distribution \21\ which galform assumes, 
together with some decrease in density away from the plane 
of the disk: 



Pi(R, z) = pi(R, 0)e 



[so Ei(i?) =2p i (i?,0)fc](24) 



Applying the local star formation assumption of \22\ i) to 
this distribution gives: 



i> a = 27r\/Gp C oid(0,0)5 / e'^dz e ^ R dR 

J — oo JO 

= \^\ m::,j 7 , . (25) 



/2 

ld \/ 47T/17-2 



Comparison between this and the usual prescription in gal- 
form exposes the big physical difference between the two 
star formation models; \22\ i) considers the gravitational col- 
lapse of the gas due to its own gravity while ( 22 a) includes 



the gravity of all component^] including the stars and dark 
matter (Af to t rather than just M co id). 

The derived analytic equivalent ( 25 I is shown in Fig[7] 
(using the fitted values of td from \ 3.4 and the illustrative 
approximatiorj^jof a constant scaleheight, h = 200pc) along- 
side the full simulation results. Discrepancies are larger at 



7 An additional factor involving the circular speed appears in 
(4.11) of Cole et al. (2000). This is omitted here (equivalent to 
setting a* =0). 

8 The dark matter does influence the creation of stars due to its 
effect on the overall stability, but this is a separate part of the 
analytic treatment (i |3.9[ l. 

9 This choice is not too critical since h appears only to the one 
half power in eqn. |25|. 
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Figure 7. The star formation efficiency from the simulation is 
shown as a solid line. The points show the result of applying the 
corrected star formation rate {25} to the disk scalelengths, rjj 
which best fit the simulation (Fig[4j. The illustrative approxi- 
mation of a constant scaleheight, h = 200pc, is used. The feint 
dotted line shows the efficiency found in empirical studies, which 
is given by eqn. {22})) with e* = 0.017 (Kennicutt 1998). 



0.1 

© 
S 

o 
O 

\ 0.05 



- 


I 

I 

I 


I 


1 




/' 

' N 
<l\ 


ih 


Hot disk gas _ 




" / \ 

I / 
l I 


i; \ 


M. t h „ t " 




I ; 
/ / 


1 


1 



5 10 

Time / Gyr 



Figure 8. A simple illustration of the correlation between star 
formation rate and the mass of he ated gas present in the disk. The 
effective timescale from equation {28}, th ea t = 7 X 10 6 years. The 
inset panel gives an indication of the almost negligible fraction of 
the total disk gas which this represents. 



early times when the assumed cold gas gas density distribu- 
tion is tenuous, but the match improves as the disk settles 



and is excellent at late times, showing that ( 25 1 is as consis- 



tent with the assumptions in Gasoline as is possible within 
this particular analytic framework. 



3. 7. 1 Burst star formation 

During merger events, and if the disk is deemed to be un- 
stable (53.91, cold disk gas in the galform model is as- 
sumed to funnel to the centre of the new system and be 
converted to stars within a timestep, which in this case is 
about 0.2Gyr (Cole et al.2000; §4.3). The incorporation of 
this non-equilibrium effect should broadly correct the total 
star formation rate on occasions when the assumed den- 



sity profile (21 1 is not such a good description of the true 



gas distribution. This is relevant here in the early stages of 
the simulation (as apparent when comparing the points and 
solid line in Fig. [7}. 



3.7.2 Schmidt- Kennicutt Law 



Equation ( 22 1) can be recognised as a version of the estab- 



lished empirical relation between star formation rate and gas 
supply (Kennicutt 1998), for which the constant of propor- 
tionality has the value e* ~ 0.017. The corresponding star 
formation efficiency for this simulated system can be cal- 
culated from its approximate half-mass radius and circular 
velocity (Fig. 10 \ and this is shown in Fig. [7] as a guide. 



3.8 Feedback 

3. 8. 1 Heating of disk gas 

The MW1 simulation treats each type II supernova event as 
a blast wave which expands to a maximum radius, r max , and 
the enclosed volume is assumed to be unable to cool for as 
long as the shell can be expected to survive, t max . If the star 
formation rate varies little over this timescale, and the shell 
subsequently cools rapidly, this will correspond to a heated 
mass, 



Mheat(t) 



rnsN 



( 47tp c 



>ld^*max 



M*(t)t h eat 



(26) 



where msN is the average mass of stars formed per super- 
nova. Now, ( 26 1 conceals the fact that both r max and £ max 



(and hence £h ea t) are, rightly, functions of the local gas den- 
sity, p and pressure, P. The precise dependence assumed is 
found in equations 9 and 11 of Stinson et al. (2006). Com- 



bining this into ( 26 1 gives 

^he a t 



/ -0.16 d -0.2n3 0.34 D -0. 7 
P(P P ) P P 



(27) 



Referring to the phase diagram in figure 4 of Stinson 
et al. (2006), gas in star forming regions appears to be ap- 
proximately isothermal at T « 10 4 K (which is expected due 
to the steep gradient in the cooling curve at this temper- 
ature). This leads to the suggestion that the environmen- 
tal dependence in (27 1 will be minimal, the pressure and 



density terms almost cancelling for constant temperature. 
Though tenuous, this provides at least a rough explana- 
tion for the extremely simple dependence which does indeed 
emerge from the full calculation, as shown in Fig. [8] 



M heat (£) w M+(t)t hea . t 



the 



7 x 10 6 yrsl (28) 



The precise definition of Afh oa t used in the figure is the mass 
of disk gas which fails to satisfy the temperature criterion for 
star formation (T < 30,000K). Thus, the integrated effect 
of feedback will be to slightly reduce the fraction of the total 
disk gas which is available for star formatiorj^] 

^ _ M co M — A-fhcat _ Af co ld ^g-j 



where tv is the particular timescale derived from equation 
(25 1. This revision is duly applied in galform, being a fair 
reflection of the physical assumptions in the simulation. 



10 For a massive galaxy such as this, the star-forming timescale is 
typically > 10 s years, so the correction to allow for re-heating 
actually makes very little difference. Quick examination of the 
inset to Fig. [8] confirms this; the fraction of the total gas which 
gets heated above the star-forming limit is tiny. 
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Figure 9. The dimesionless ratio which is used as a stability 
critereon in the galform model, shown for this particular re- 
alisation (solid line) and for a relisation where the collapse of 
"unstable" disks has not been enforced (dotted line). 



3.8.2 Outflow 

To discover the extent of outflow from the galaxy in the 
simulation, all particles were identified which have moved 
from the inside to the outside of a radius of 30 comoving 
kpc. The total mass of such material turned out to be even 
less than the heated disk mass (see Fig. [8| and can safely be 
ignored for the purposes of this study. The outflow of gas as 
a result of supernovae in the central galaxy is therefor set 
to zero in this version of galform, for consistency with the 
Gasoline simulation. 

It is important to note that this feedback model is more 
effective at lower halo masses, due to shallower potential 
wells (Brooks et al. 2007). This dependence is also accounted 
for in the usual galform feedback model, where the rate 
of mass outflow is assumed to be at least equal to the star 
formation rate, and thousands of times higher for smaller 
systems (Bower et al. 2006). The consequence of this differ- 
ence is investigated in §4.3| 



3.9 Disk stability 

The analytic condition which is used to assess the stability 
of disks in the galform model follows Efstathiou, Lake & 
Negroponte (1982). It deems disks to be unstable to bar 
formation if: 



GMdisk/rD 



> 1.56 



(30) 



When this inequality is met, all disk material is transferred 
to a central bulge. 

In the usual galform model, some set fraction of disk 
mass is also transferred to a central black hole (Bower et al. 
2006). This was used to arrive at the value in (301 which 



was found to best match the Magorrian relation between 
bulge mass and black hole mass, A^bh ~ Af^ u 5g C , observed 
by Haring & Rix (2004). In this study, the black hole ac- 
cretion fraction is set to zero for the purposes of consitency 
with the simulation, where this aspect of evolution is not 
considered. Being physically well motivated, the condition 
( |30[ ) has been applied throughout this study. Conveniently, 
the precise formalism used need not cause too much concern 
here; the change that its application produces in the compo- 
nent masses (in Fig. |11| for example) is barely enough to be 
noticed. This can be understood from Fig. [5] which shows 
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Figure 10. The evolution of the disk rotation velocity, v c , (top 
panel) and the disk scalelength, rrj, (lower panel) for both realisa- 
tions. The circular velocity of the simulated disk (black squares) 
is taken at 5 comoving kpc from the centre. The approximate gas 
scalelengths (blue squares) are found by fitting t ]2 1 [ ) to the simu- 
lated gas distributions in Fig. |4| T he stellar scalelengths (red tri- 
angles) are taken from fitting pi) the the stellar profiles, shown 
later in Fig. |B1| Dashed lines are the values calculated in gal- 
form from conservation of angular momentum (j |3.6| - gas and 
stars are assumed to have the same scalelength). 



that only a very limited adjustment to the disk mass is re- 



quired to prevent the limit in ( 30 1 from being reached. The 



mass aggregation of this particular system, and the fact that 
the specific angular momentum of the gas is relatively higlj^] 
(A > 0.07 for the latter half of the halo's history) leads to 
the stabilisation of the system, even if the redistribution of 
mass is not enforced. 



4 RESULTS 

4.1 Matched Model 

4-1.1 Structure 

The size and rotation speed are shown for both realisations 
in Fig. |10| The predicted galform value is based on con- 
serving the initial angular momentum ( |19[ ) of in-falling gas; 
an assumption which is particularly effective in this case 
where the distributions for the initial and final mass were 
appropriate (Fig|4|. This correspondence is encouraging for 
the approach of calculating disk sizes from the spin param- 
eter, A, of the halo from which the gas cooled. 



11 The criterion l |30| can also be written: 
Af disk > 1.2[A/ DM (r 1/2 ) + M sph (r 1/2 )] , 



(31) 



where Mdm and M sp ^ are the masses of the dark matter and 
spheroid components. The disk's eventual half mass radius will 
enclose more dark matter if the original spin parameter is higher. 
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Figure 11. The evolution of hot halo gas, cold disk gas and stars 
according to both modeling approaches. The star formation rate 
in galform is calculated using assumptions identical to those 
in the simulation l |22^ t), but taking into account their integrated 
effect over the whole disk profile 1211. 



4.1.2 Cooled Mass 

The total cooled mass (stars and cold gas), calculated using 
the methods of §3.2| can be seen from Fig. |11| Agreement 
with the simulation is excellent, given the inherent difference 
in the methods of calculation, and bolsters confidence in the 
rather intricate calculation of the energy radiated by the 
system throughout its complex merger history 1 13)18 1. 

The division of the cooled mass into stars and gas are 
also shown in Fig. |11| The evident agreement can be under- 
stood by reference to Fig's |4j [7] and [8] which demonstrate 
that the considerable complexities of the simulation reduce 
to relatively simple relationships when integrated over the 
entire disk. 



4-1.3 Discrepancies 

There is a minor discrepancy at early times (t < 8Gyr) 
in the cooled mass generated by the two modelling tech- 
niques. This is unlikely to be due to differences in cooling 
from the hot halo: given the similar density distributions 
(Fig. |4| and cooling function, the two models should yield 
similar amounts of cooled gas, all else being equal. 

The discrepancy is also not due to a difference in the 
time of development of a hot gas halo. The top panel of Fig. 
4 shows that accreted gas in the galform model transitions 
from unshocked, cold gas to shocked, hot gas at an age of 
about 2 Gyr, while the bottom panel shows that this time is 
also when the gasoline model begins to develop a hot gas 
halo. 

This demonstrates that the adopted analytic separation 
of unshocked and shocked gas (free fall limited regime versus 
cooling limited regime) is an excellent analog to the transi- 
tion between cold gas accretion and shocked gas accretion 
seen in simulations as a galaxy crosses the threshold mass ca- 
pable of developing a hot halo. It is important to note that 



10 



Standard GASOLINE 
Forced shock-heatine 




10 



Time / Gyr 



Figure 12. Having focussed on modifications to galform, this 
figure shows the effects of modifying the SPH code gasoline to 
emulate the semi-analytic technique. The solid black line is the 
star formation rate (lower panel) and its integral (upper panel) 
which occurs in the standard gasoline simulation which has been 
under discussion. The dotted red lines show the consequence of 
delaying star formation from cold accreted gas is delayed by a 
cooling time. 



this transition has always been included in analytic mod- 
els of galaxy formation, and thus cold gas accretion at low 
galaxy mass does not alter the standard picture of galaxy 
formation. 

The discrepancy in cooled mass arises after the devel- 
opment of this hot halo. The simulated galaxy continues to 
accrete cold gas via filaments that penetrate within the hot 
gas halo (see Fig. [3] here, figure 5 of Brooks et al. (2009) 
and also Ocvirk et al. (2008), Agertz et al. (2009) and Dekel 
et al. (2009)), while cold gas accretion ends in the analytic 
model. 



4.2 Adapting gasoline to emulate semi-analytics 

Brooks et al. (2009) demonstrated that the inclusion of cold 
gas accretion along filaments leads to an earlier phase of star 
formation than predicted if all gas was initially shock heated, 
due to the shorter cooling times onto the central galaxy of 
the cold gas. Fig. 12 shows that the simulated galaxy does 
indeed form stars earlier than the analytic model, though 
the overall discrepancy is never more than a factor of two. 

In an attempt to adapt the simulations to include simi- 
lar physics as the analytic model, Fig. |12| shows the effect on 
the stellar growth in the simulations if all of the cold gas had 
instead been shocked. The hot gas has cooling times several 
Gyr longer than the cold gas. This is quantified, and a delay 
has been added to the formation time of the stars spawned 
from the cold accreted gas (Brooks et al. 2009). As seen in 
Fig. |12[ this delay leads to a slower build of up stellar mass 
in the simulations. Comparison with Fig]ll| shows this to 
be in even better agreement with the stellar growth of the 
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Figure 13. The results of applying the unadjusted version of 
galform (Bower et al. 2006) to the same merger tree as used 
for Fig. |ll| and indeed the whole paper. To bring out the contrast 
between the two models, the shading now represents the galform 
predictions and the lines show the simulation results. The upper 
panel is the equivalent to Fig. [9] showing the stability of the disk 
in this version for this model. 



analytic model, which neglects this cold gas accretion along 
filaments. 



4.3 Comparison with standard galform 
assumptions 

In order to represent the same physical assumptions as this 
particular simulation, the semi-analytic code had to be sig- 
nificantly modified from the version which has been used 
most recently to study the collective properties of galaxy 
samples (Bower et al. 2006; Stringer et al. 2009). There 
are three very significant differences between this version of 
galform and the simulation studied here, all of which are 
manifest in the comparison of the history of the simulation 



and the Bower et al. (2006) model (Fig. 13 1 



4-3.1 Inclusion of AGN feedback 



following criteria are satisfied, it is assumed that no hot halo 
gas will be able to cool onto the disk. 

(1) The resulting Eddington luminosity of the black hole is 
at least 4% of the cooling luminosity of the halo. 

(2) The cooling time Q exceeds the freefall time |6| by 
some chosen factor (in this case, t coo \ > 0.58tfj). 

This effect is visible in Fig.[l3]as a period of rapidly decreas- 
ing cold gas mass; during such phases it is no longer being 
replenished by gas cooling in from the halo. 



4-3.2 Stronger supernova feedback 

This is assumed to be extremely strong in the Bower et al. 
(2006) model, the justification being that a more modest 
conversion of supernova energy to gas outflow would allow 
the formation of too many low-mass galaxies. (Indeed, it 
appears that the predicted number may still be too high 
(Stringer et al. 2009) even with this generous allocation of 
feedback energy.) 

In conjunction with the inclusion of AGN heating, this 
assumption leads to a system with a much greater hot gas 
component: As soon as the supply of gas to the disk can no 
longer be replenished by cooling, all the gas in the system 
quickly ends up in the hot component. This effect is ad- 
ditionally enforced by the low star formation efficiency (see 
below) , which prevents disk mass being locked up into stars. 



4-3.3 Inefficient star formation 

The star formation assumptions in Gasoline and the Cole 



et al. (2000) model have already been contrasted in 53.7 



When there are instabilities (£3.9) or merger events in this 
model which trigger disk collapse, 5% of the available disk 
gas is assumed to accrete onto a central black hole. If the 



The assumed efficiency in the Bower et al. (2006) model is 
to be so low (e* = 0.0029) that, even after allowing for the 
structure and gas fraction, the final conversion rate is about 
a factor of 10 lower than in the simulation. (M*/M co id ~ 
to 0.2/Gyr as opposed to the rate of 0.5 to 3/Gyr seen in 
Fig.0. 

This stands out in Fig |13| The large mass of cold gas 
which is supported by the disk in the semi-analytic realisa- 
tion (except when cooling is shut off by AGN) cannot be 
sustained in the simulation due to much more effective con- 
version to stars. 



5 SUMMARY 

The merger history of a simulated, Milky Way-like galaxy 
has been used to recreate the same system using the gal- 
form semi-analytic model. 

The incoming trajectories of satellite halos were inves- 
tigated and found to be consistent with the usual, assumed 
distributions. Information on the different modes of accre- 
tion onto the main halo is shown to be inherently contained 
in the merger tree (and would be similarly contained in a 
statistically generated merger tree which may be used in the 
study of larger volumes) . 

The subsequent fate of accreted gas is not modelled 
equivalently by the two methods, but the analytic framework 
does account for both shocked and unshocked gas. Further- 
more, the estimated formation time of the hot halo is in 
excellent agreement with the simulation. 
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The single-parameter analytic distributions which are 
assumed in galform are found to be a satisfactory descrip- 
tion of the radial distribution of gas and stars inside the 
simulated halo. It was also shown that the choice of hot 
halo scalelength is of minor importance in the calculation of 
cooled gas mass. 

However, the enforcement of spherical symmetry in 
these analytic distributions means that any subsequent cold 
accretion along filaments (persisting in spite of shocks else- 
where) is neglected. Thus, the total quantity of unshocked 
gas onto the central galaxy is underestimated, as expected 
by Brooks et al. (2009). Some gas is consequently delayed in 
its arrival at the central galaxy, though it is found to proceed 
at later times. 

The disk structure is predicted analytically by conserv- 
ing the angular momentum of the cooled gas. Resulting 
scalelengths are within a factor of two of the simulation re- 
sults and the circular velocities even closer. This is quite 
satisfactory for the case of a single halo, but more work is 
needed to establish whether the two techniques will con- 
sistently agree to this level for an arbitrary set of initial 
conditions (and whether any disagreement is biased). 

With regard to the subsequent evolution of the disk, 
the galform code was altered to apply the same funda- 
mental assumptions as the simulation by considering their 
integrated effect over the assumed disk structure. Subsequent 
agreement was excellent given this immense reduction in 
complexity from many thousands of SPH particles to a few 
one-dimensional equations. 

So, by assuming the same physics and using the same 
initial conditions, the two techniques predict a final system 
which is recognisably the same, despite fundamental differ- 
ences in the way that these assumptions are applied and the 
evolution followed. This suggests that equations which at- 
tempt to characterise the emergent behaviour of the system 
may, with sufficient understanding, become as reliable as 
those which emulate the properties of particles themselves. 

With such potential consistency thus established, the 
same system was then evolved under physical assumptions 
which have been more frequently adopted within galform 
(Bower et al. 2006). The resulting system (at any time) is 
not recognisably the same as that predicted by the simula- 
tion: having less than a quarter of the stellar mass, barely 
half the rotational velocity and with negligible star forma- 
tion at z=0. Though not all aspects of a system need to be 
correctly represented to learn from the modelling exercise, 
this divergence from the same initial conditions proves that 
at least one of these published models poorly represents the 
true nature of galaxy formation. 

It is hoped that this detailed comparison between two 
very different theoretical approaches will promote the aware- 
ness of their respective limitations and help lead them both 
towards a truer description of real galaxies. 
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APPENDIX A: SATELLITES 

The orbits of sub halos around their host can be charac- 
terised (Lacey & Cole 1993) by the dimensionless parameter 

e=|^!°' 7 Y^Y , (Ai) 



ME) 



where r c (E) and j c (E) are the radius and specific angular 
momentum of a circular orbit with energy E. For a potential 
with a flat rotation curve; where the circular velocity at all 
radii is ti c , 



$(r) 



$(r v ) +v 2 c In ( — 



(A2) 



Hence, for a subhalo at position, r with velocity v, the radius 
of a circular orbit with the same energy is: 



r c (v, r) = | r | exp 



v 

2vl 



(A3) 



The orbital parameters are found by simply applying equa- 



tions ( Al I to ( A3 1 to the position and velocity vectors of 



each satellite halo (after moving to a frame where these vec- 
tors are both for the host halo): 

(, n 70 / ' * — 1 ^ \ 1-22 
v c r v J 



ME) 
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Figure Al. Left panel: The velocities and radii of the satellite 
halos in the simulation of Brooks et al. (2009), at the timestep 
immediately before they are deemed to have merged with their 
host. The dots are those which appear in the distribution of 
in Fig. [I] crosses are those that have > 4.5 and hence do not 
appear. The line shows the constant value of © corresponding 
to that cut-off, for the given approach angle. Right panel: The 
distribution of approach angles of these satellites alongside the 
distribution found from the study by Benson (2005). 
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The galform model assigns each satellite a value for this 
parameter, Q, drawn at random from the log-normal distri- 
bution that it was found to follow in simulations by Tormen 
(1997): 

dlnn 



dln9 



exp 



(lne + 0.14) 2 
2^2 



(A5) 



with a = 0.26. The sample of orbit parameters produced by 
this random assignment is shown in Fig. |A1| together with 
the actual values from the Brooks et al. (2009) simulation. 

The fraction of the halo free-fall time that it takes for 
the satellite to merge through dynamical friction is assumed 
to be proportional to O: 



[G Pv (z)Y 



e 



0.3722 



ln(Mhost/M sat ) M si 



(A6) 



This is based on the standard Chandrasekhar formula for the 
dynamical friction and appears originally in Lacey & Cole 
(1993). p v is the mean density at the virial radius (deter- 
mined by the cosmology, not the specific halo's properties). 

The distribution of angles at which sub structures en- 
ter their host halo has been studied by Benson (2005) for 
the VLS and VIRGO simulations and and found to have 
a repeatable distribution. This distribution can be applied 
in the galform model to generate an alternative set of or- 
bital parameters, G, which are shown in Fig. |A1| alongside 
the standard assumption of the log-normal distribution, and 
the results from gasoline. 

Merger times predicted by ( |A6[ ) have been compared by 
Jiang et al. (2008) with the results of GADGET2 simulations 
(Springel 2005). The agreement found was of the order of a 
factor of two, which is deemed acceptable for the continued 
use of this formula in galform. However, since (A6l was 



found to consistently underestimate the simulated merger 
time, the improved fitting formula proposed byJiang et al. 
(2008) may well be adopted in future. 
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Unfortunately, due to the very different definition of 



merging adopted by the halo finder used here (63.1 1 , a mean- 



ingful comparison of the respective times in the two realisa- 
tion has not been possible. 



APPENDIX B: DISK PROFILES 

The stellar disk radii that appear in Fig. [5] are generated 
from analysis of the distribution of stars in the simulation. 
Fig. |B1| compares the distribution of stellar mass from the 
simulation with the analytic forms assumed by galform, 



both for disk stars (21 1 and for the mass of stars in the 



bulge, assumed to be distributed such that the projected 
surface density profile is given by: 



^buige oc exp 



(Bl) 



Though the agreement is not comprehensive, the ana- 
lytic forms are an adequate description of the global distri- 
bution of simulated stars for the majority of the systems evo- 
lution. Unsurprisingly, the simple profiles of (21 1 and ( Bl I 



fail to describe the simulated system at early times, before 
an ordered galactic system has formed. 

The fact that galform adheres, at all times, to the fit- 
ting forms that are really only appropriate to recent times 
(z < 2) is clearly a valid criticism of this technique. How- 
ever, as long as such forms are indeed a good description 
of low-redshift systems, then the only way their earlier mis- 
judgements can significantly mislead the predictions of the 
model is in the early time star formation rates (Fig[7|, which 
are based on this assumed structure. 




(l + z) r / kpc 

Figure Bl. The stellar mass profile from the simulation of 
Brooks et al. (2009) at six different redshifts (solid black lines). 
The dotted blue lines show the analytic forms which are assumed 
by the galform model for the disk f21\ . Dotted red lines show 
the bulge components flBlfl . The total stellar contribution, shown 
by dotted magenta lines, is constrained to integrate to the same 
mass as the soild lines. The two lines are therefore fitted using 
three free parameters: The relative fraction of mass in the two 
components, the disk scalelength, rjy and the scalelength of the 
bulge profile, r^. 



